/*******************************************************************************
	   Authors: Kaixin Wang
	   Project: ETS
	   Purpose: Plot time series of CEMS data availability and PM emissions

	   Date created: 18 February 2022
	   Version:      STATA 17 MP
	   
	   Last edited: 2 June 2022
	   Edited by: 	Kaixin		
   ****************************************************************************/

set more off 
clear matrix
clear all
pause on

********************************************************************************
************************** [1] PLOT DATA AVAILABILITY **************************
********************************************************************************

use "$IMPUTATION_DATA_OUT/PlantWeekPMMassNoRule.dta", clear
local n_ctrl = 136
local n_trt = 156
local n_tot = 292

gen D_report = 0
replace D_report = 1 if ind_week_mass_val_rule0 != .
bysort week treatmentstatus: egen weekly_pct_report = mean(100 * D_report)
bysort week treatmentstatus: egen weekly_mean_PM = mean(ind_week_mass_val_rule0)
keep week treatmentstatus weekly_pct_report weekly_mean_PM
duplicates drop
gen week_new = td(31dec2018) + (week - 1) * 7
format week_new %tddd-Mon-YY
gen week_plot = week_new + 3
format week_plot %tddd-Mon-YY

local xtext_pos = 0
local x_opt xtick(21655 21685 21716 21746 21777 21808 21838 21869 21899 21930 21961 21990 22021 22051 22082 22112 22143 22174 22204 22235 22265 22296 22327 22355 22386) ///
	xlabel(21655 21746 21838 21930 22021 22112 22204 22296 22386, angle(35) labsize(medsmall)) ///
	xline(21745, lwidth(vthin) lcolor(black) lpattern(tight_dot)) ///
	text(`xtext_pos' 21746 "Mock-I", place(ne) orientation(vertical) size(small) color(gs8)) ///
	xline(21774, lwidth(vthin) lcolor(black) lpattern(tight_dot)) ///
	text(`xtext_pos' 21775 "Mock-II", place(ne) orientation(vertical) size(small) color(gs8)) ///
	xline(21808, lwidth(vthin) lcolor(black) lpattern(tight_dot)) ///
	text(`xtext_pos' 21809 "Period-I", place(ne) orientation(vertical) size(small) color(gs8)) ///
	xline(21823, lwidth(8.2) lcolor(gs14%30) lpattern(solid)) ///
	xline(21838, lwidth(vthin) lcolor(black) lpattern(tight_dot)) ///
	text(`xtext_pos' 21839 "Period-II", place(ne) orientation(vertical) size(small) color(gs8)) ///
	xline(21853.5, lwidth(8) lcolor(gs14%30) lpattern(solid)) ///
	xline(21869, lwidth(vthin) lcolor(black) lpattern(tight_dot)) ///
	text(`xtext_pos' 21870 "Period-III", place(ne) orientation(vertical) size(small) color(gs8)) ///
	xline(21892, lwidth(12.5) lcolor(gs14%30) lpattern(solid)) ///
	xline(21915, lwidth(vthin) lcolor(black) lpattern(tight_dot)) ///
	text(`xtext_pos' 21916 "Period-IV", place(ne) orientation(vertical) size(small) color(gs8)) ///
	xline(21930.5, lwidth(8) lcolor(gs14%30) lpattern(solid)) ///
	xline(21946, lwidth(vthin) lcolor(black) lpattern(tight_dot)) ///
	text(`xtext_pos' 21947 "Period-V", place(ne) orientation(vertical) size(small) color(gs8)) ///
	xline(21960.5, lwidth(7.8) lcolor(gs14%30) lpattern(solid)) ///
	xline(21975, lwidth(vthin) lcolor(black) lpattern(tight_dot)) ///
	text(`xtext_pos' 21976 "Period-VI", place(ne) orientation(vertical) size(small) color(gs8)) ///
	xline(21985.5, lwidth(5.5) lcolor(gs14%30) lpattern(solid)) ///
	xline(21996, lwidth(vthin) lcolor(black) lpattern(tight_dot)) ///
	text(`xtext_pos' 21997 "Interregnum", place(ne) orientation(vertical) size(small) color(edkblue%70)) ///
	xline(22098, lwidth(55.5) lcolor(ebg%50) lpattern(solid)) ///
	xline(22200, lwidth(vthin) lcolor(black) lpattern(tight_dot)) ///
	text(`xtext_pos' 22201 "Mock-III", place(ne) orientation(vertical) size(small) color(gs8)) ///
	xline(22231, lwidth(vthin) lcolor(black) lpattern(tight_dot)) ///
	text(`xtext_pos' 22232 "Interregnum", place(ne) orientation(vertical) size(small) color(edkblue%70)) ///
	xline(22240.5, lwidth(5) lcolor(ebg%50) lpattern(solid)) ///
	xline(22250, lwidth(vthin) lcolor(black) lpattern(tight_dot)) ///
	text(`xtext_pos' 22251 "Period-VII", place(ne) orientation(vertical) size(small) color(gs8)) ///
	xline(22265.5, lwidth(8) lcolor(gs14%30) lpattern(solid)) ///
	xline(22281, lwidth(vthin) lcolor(black) lpattern(tight_dot)) ///
	text(`xtext_pos' 22282 "Period-VIII", place(ne) orientation(vertical) size(small) color(gs8)) ///
	xline(22296.5, lwidth(8) lcolor(gs14%30) lpattern(solid)) ///
	xline(22312, lwidth(vthin) lcolor(black) lpattern(tight_dot)) ///
	text(`xtext_pos' 22313 "Period-IX", place(ne) orientation(vertical) size(small) color(gs8)) ///
	xline(22326, lwidth(7.5) lcolor(gs14%30) lpattern(solid)) ///
	xline(22340, lwidth(vthin) lcolor(black) lpattern(tight_dot)) ///
	text(`xtext_pos' 22341 "Period-X", place(ne) orientation(vertical) size(small) color(gs8)) ///
	xline(22355.5, lwidth(8) lcolor(gs14%30) lpattern(solid)) ///
	xline(22371, lwidth(vthin) lcolor(black) lpattern(tight_dot)) 	

twoway ///
	connect weekly_pct_report week_plot if treatmentstatus == "Treatment", ///
	clwidth(vthin) clcolor(midblue) mcolor(midblue) msymbol(T) msize(vsmall) || ///
	connect weekly_pct_report week_plot if treatmentstatus == "Control", ///
	clwidth(vthin) clcolor(gs6) mcolor(gs6) msize(vsmall) lpattern(shortdash) ///
	plotregion(color(white)) bgcolor(white) graphregion(color(white)) ///
	xtitle("") ///
	ytitle("Plants reporting (%)", size(medsmall) axis(1)) ///
	ylabel(0(20)100, labsize(medsmall) angle(horizontal) axis(1)) ///
	`x_opt' ///
	legend(cols(1) size(medsmall) region(lstyle(none)) position(6)) ///
	legend(on order(1 "Treatment (n=`n_trt')" 2 "Control (n=`n_ctrl')")) ///
	xsize(9) ysize(4)

graph export "$EMISSIONS_FIGS/Figure_C2.pdf", replace


********************************************************************************
************************** [2] PLOT WEEKLY EMISSION ****************************
********************************************************************************

local RuleNo 0 A B

foreach r of local RuleNo {
	
	use "$IMPUTATION_DATA_OUT/PlantWeekPMMassRule`r'.dta", clear

	gen D_report = 0
	replace D_report = 1 if ind_week_mass_val_rule0 != .
	bysort week treatmentstatus: egen weekly_pct_report = mean(100 * D_report)
	bysort week treatmentstatus: egen weekly_mean_PM = mean(ind_week_mass_val_rule0)
	keep week treatmentstatus weekly_pct_report weekly_mean_PM
	duplicates drop

	gen week_new = td(31dec2018) + (week - 1) * 7
	format week_new %tddd-Mon-YY
	gen week_plot = week_new + 3
	format week_plot %tddd-Mon-YY

	gen weekly_mean_PM_prorated = weekly_mean_PM / 7 * 30

	** Plot WEEKLY Mean PM Mass by Treatment Status (NO IMPUTATION)

	local xtext_pos = 0
	local x_opt xtick(21655 21685 21716 21746 21777 21808 21838 21869 21899 21930 21961 21990 22021 22051 22082 22112 22143 22174 22204 22235 22265 22296 22327 22355 22386) ///
		xlabel(21655 21746 21838 21930 22021 22112 22204 22296 22386, angle(35) labsize(medsmall)) ///
		xline(21745, lwidth(vthin) lcolor(black) lpattern(tight_dot)) ///
		text(`xtext_pos' 21746 "Mock-I", place(ne) orientation(vertical) size(small) color(gs8)) ///
		xline(21774, lwidth(vthin) lcolor(black) lpattern(tight_dot)) ///
		text(`xtext_pos' 21775 "Mock-II", place(ne) orientation(vertical) size(small) color(gs8)) ///
		xline(21808, lwidth(vthin) lcolor(black) lpattern(tight_dot)) ///
		text(`xtext_pos' 21809 "Period-I", place(ne) orientation(vertical) size(small) color(gs8)) ///
		xline(21823, lwidth(8) lcolor(gs14%30) lpattern(solid)) ///
		xline(21838, lwidth(vthin) lcolor(black) lpattern(tight_dot)) ///
		text(`xtext_pos' 21839 "Period-II", place(ne) orientation(vertical) size(small) color(gs8)) ///
		xline(21853.5, lwidth(8) lcolor(gs14%30) lpattern(solid)) ///
		xline(21869, lwidth(vthin) lcolor(black) lpattern(tight_dot)) ///
		text(`xtext_pos' 21870 "Period-III", place(ne) orientation(vertical) size(small) color(gs8)) ///
		xline(21892, lwidth(12) lcolor(gs14%30) lpattern(solid)) ///
		xline(21915, lwidth(vthin) lcolor(black) lpattern(tight_dot)) ///
		text(`xtext_pos' 21916 "Period-IV", place(ne) orientation(vertical) size(small) color(gs8)) ///
		xline(21930.5, lwidth(8) lcolor(gs14%30) lpattern(solid)) ///
		xline(21946, lwidth(vthin) lcolor(black) lpattern(tight_dot)) ///
		text(`xtext_pos' 21947 "Period-V", place(ne) orientation(vertical) size(small) color(gs8)) ///
		xline(21960.5, lwidth(7.8) lcolor(gs14%30) lpattern(solid)) ///
		xline(21975, lwidth(vthin) lcolor(black) lpattern(tight_dot)) ///
		text(`xtext_pos' 21976 "Period-VI", place(ne) orientation(vertical) size(small) color(gs8)) ///
		xline(21985.5, lwidth(5.4) lcolor(gs14%30) lpattern(solid)) ///
		xline(21996, lwidth(vthin) lcolor(black) lpattern(tight_dot)) ///
		text(`xtext_pos' 21997 "Interregnum", place(ne) orientation(vertical) size(small) color(edkblue%70)) ///
		xline(22098, lwidth(55) lcolor(ebg%50) lpattern(solid)) ///
		xline(22200, lwidth(vthin) lcolor(black) lpattern(tight_dot)) ///
		text(`xtext_pos' 22201 "Mock-III", place(ne) orientation(vertical) size(small) color(gs8)) ///
		xline(22231, lwidth(vthin) lcolor(black) lpattern(tight_dot)) ///
		text(`xtext_pos' 22232 "Interregnum", place(ne) orientation(vertical) size(small) color(edkblue%70)) ///
		xline(22240.5, lwidth(5) lcolor(ebg%50) lpattern(solid)) ///
		xline(22250, lwidth(vthin) lcolor(black) lpattern(tight_dot)) ///
		text(`xtext_pos' 22251 "Period-VII", place(ne) orientation(vertical) size(small) color(gs8)) ///
		xline(22265.5, lwidth(8) lcolor(gs14%30) lpattern(solid)) ///
		xline(22281, lwidth(vthin) lcolor(black) lpattern(tight_dot)) ///
		text(`xtext_pos' 22282 "Period-VIII", place(ne) orientation(vertical) size(small) color(gs8)) ///
		xline(22296.5, lwidth(8) lcolor(gs14%30) lpattern(solid)) ///
		xline(22312, lwidth(vthin) lcolor(black) lpattern(tight_dot)) ///
		text(`xtext_pos' 22313 "Period-IX", place(ne) orientation(vertical) size(small) color(gs8)) ///
		xline(22326, lwidth(7) lcolor(gs14%30) lpattern(solid)) ///
		xline(22340, lwidth(vthin) lcolor(black) lpattern(tight_dot)) ///
		text(`xtext_pos' 22341 "Period-X", place(ne) orientation(vertical) size(small) color(gs8)) ///
		xline(22355.5, lwidth(8) lcolor(gs14%30) lpattern(solid)) ///
		xline(22371, lwidth(vthin) lcolor(black) lpattern(tight_dot)) 

	twoway ///
		connect weekly_mean_PM_prorated week_plot if treatmentstatus == "Treatment" & week_new < date("2020mar23", "YMD"), ///
		clwidth(vthin) clcolor(midblue) mcolor(midblue) msymbol(T) msize(vsmall)|| ///
		connect weekly_mean_PM_prorated week_plot if treatmentstatus == "Treatment" & week_new > date("2020oct5", "YMD") & week_new < date("2020nov9", "YMD"), ///
		clwidth(vthin) clcolor(midblue) mcolor(midblue) msymbol(T) msize(vsmall) || ///
		connect weekly_mean_PM_prorated week_plot if treatmentstatus == "Treatment" & week_new > date("2020nov30", "YMD"), ///
		clwidth(vthin) clcolor(midblue) mcolor(midblue) msymbol(T) msize(vsmall) || ///
		connect weekly_mean_PM_prorated week_plot if treatmentstatus == "Control" & week_new < date("2020mar23", "YMD"), ///
		clwidth(vthin) clcolor(gs6) mcolor(gs6) msymbol(o) msize(vsmall) lpattern(shortdash) || ///
		connect weekly_mean_PM_prorated week_plot if treatmentstatus == "Control" & week_new > date("2020oct5", "YMD") & week_new < date("2020nov9", "YMD"), ///
		clwidth(vthin) clcolor(gs6) mcolor(gs6) msymbol(o) msize(vsmall) lpattern(shortdash) || ///
		connect weekly_mean_PM_prorated week_plot if treatmentstatus == "Control" & week_new > date("2020nov30", "YMD"), ///
		clwidth(vthin) clcolor(gs6) mcolor(gs6) msymbol(o) msize(vsmall) lpattern(shortdash) || ///
		pci 1728 21745 1728 21837, lcolor(red) lwidth(vthin) || ///
		pci 1235 21837 1235 21868, lcolor(red) lwidth(vthin) || ///
		pci 1111 21868 1111 21914, lcolor(red) lwidth(vthin) || ///
		pci 1049 21914 1049 21995, lcolor(red) lwidth(vthin) || ///
		pci 1049 22200 1049 22230, lcolor(red) lwidth(vthin) || ///
		pci 1049 22250 1049 22371, lcolor(red) lwidth(vthin) ///
		plotregion(color(white)) bgcolor(white) graphregion(color(white)) ///
		xtitle("") ///
		ytitle("Mean PM mass (kg/month)", size(medsmall) axis(1)) ///
		ylabel(0(500)3000, labsize(medsmall) angle(horizontal)) ///
		`x_opt' ///
		text(1658 21745 "1728 kg", place(ne) size(tiny) color(red)) ///
		text(1265 21837 "1235 kg", place(ne) size(tiny) color(red)) ///
		text(1041 21868 "1111 kg", place(ne) size(tiny) color(red)) ///
		text(979 21914 "1049 kg", place(ne) size(tiny) color(red)) ///
		text(979 22200 "1049 kg", place(ne) size(tiny) color(red)) ///
		text(979 22250 "1049 kg", place(ne) size(tiny) color(red)) ///
		legend(cols(1) size(medsmall) region(lstyle(none)) position(6)) ///
		legend(on order(4 "Control (n=`n_ctrl')" 1 "Treatment (n=`n_trt')")) ///
		xsize(9) ysize(4) 

		if "`r'" == "0" local figure_name = "Figure_5.pdf" 
		if "`r'" == "A" local figure_name = "Figure_C3_A.pdf" 
		if "`r'" == "B" local figure_name = "Figure_C3_B.pdf" 
		
	graph export "$EMISSIONS_FIGS/`figure_name'", replace
}
